Petri net model using AlgebraicPetri.jl

Micah Halter (@mehalter), 2020-07-13

Introduction

One representation of the SIR model is to think of it as the combination of two interactions, transmission and recovery. AlgebraicPetri.jl allows you to define a theory of modeling (e.g. Epidemiology), and then provides a DSL for defining models in that theory as open dynamical systems where the underlying theory is preserved. This implementation defines the SIR model as the composition of two interactions defined at domain-level semantics, transmission and recovery, and then generates an appropriate Petri.jl model that can generate ODE, SDE, and jump process solvers.

Libraries

using AlgebraicPetri.Epidemiology
using Petri
using Catlab.Theories
using Catlab.CategoricalAlgebra.ShapeDiagrams
using Catlab.Graphics
using OrdinaryDiffEq
using StochasticDiffEq
using DiffEqJump
using Random
using Plots

# helper function to visualize categorical representation
display_wd(ex) = to_graphviz(ex, orientation=LeftToRight, labels=true);
display_wd (generic function with 1 method)

Define a Theory

AlgebraicPetri comes packaged with an Epidemiology module with a basic, predefined theory of epidemiology models. The source starts by defining a theory in a specified category. Here we have 4 types ($S$: susceptible, $E$: exposed, $I$: infected, $R$: recovered, $D$: dead) and 5 domain processes ($transmission: S \otimes I \rightarrow I$, $exposure: S \otimes I \rightarrow E \otimes I$, $illness: E \rightarrow I$, $recovery: I \rightarrow R$, $death: I \rightarrow D$). This defines a theory thata can be use to describe epidemiology models at domain-level semantics.

@present InfectiousDiseases(FreeBiproductCategory) begin
    S::Ob
    E::Ob
    I::Ob
    R::Ob
    D::Ob
    transmission::Hom(SI, I)
    exposure::Hom(SI, EI)
    illness::Hom(E,I)
    recovery::Hom(I,R)
    death::Hom(I,D)
end

To be able to build up large Petri nets using these terms, we must define the building block Petri nets that describe each of the interactions. For this we have the following code that defines 3 petri nets: spontaneous_petri defines a petri net with 2 states, and a transition from $S_1$ to $S_2$ to represent a spontaneous change; transmission_petri defines a petri net where $S_1$ and $S_2$ interact and produce 2 $S_2$; and lastly exposure_petri where there are 3 states and a transition where $S_1$ and $S_2$ interact and produce an $S_2$ and $S_3$ to represent an infected population causing the susceptible population to become exposed. Lastly, we need to define a dictionary that connects the objects from the theory to the petri net objects that define the interactions.

ob = PetriCospanOb(1)
spontaneous_petri = PetriCospan([1], Petri.Model(1:2, [(Dict(1=>1), Dict(2=>1))]), [2])
transmission_petri = PetriCospan([1], Petri.Model(1:2, [(Dict(1=>1, 2=>1), Dict(2=>2))]), [2])
exposure_petri = PetriCospan([1, 2], Petri.Model(1:3, [(Dict(1=>1, 2=>1), Dict(3=>1, 2=>1))]), [3, 2])

const FunctorGenerators = Dict(S=>ob, E=>ob, I=>ob, R=>ob, D=>ob,
        transmission=>transmission_petri, exposure=>exposure_petri,
        illness=>spontaneous_petri, recovery=>spontaneous_petri, death=>spontaneous_petri)

Transitions

Using the categorical framework provided by the AlgebraicJulia environment, we can think of building models as the combination of known building block open models. For example we have $transmission: S \otimes I \rightarrow I$ and $recovery: I \rightarrow R$ interactions defined in the Epidemiology module of AlgebraicPetri which can be visualized as the following Petri nets.

Transmission (where $S_1 = S$ and $S_2 = I$):

Graph(decoration(F_epi(transmission)))

Recovery (where $S_1 = I$ and $S_2 = R$):

Graph(decoration(F_epi(recovery)))

In this approach we can think of an sir model as the composition of transmission and recovery. This allows us to define the relationship that the infected population coming out of the transmission interaction is the same as population of infected in the recovery interaction.

sir_wiring_diagram = transmission  recovery
display_wd(sir_wiring_diagram)

using the function F_epi provided by AlgebraicPetri.Epidemiology, we can convert this categorical definition of SIR to the Petri net representation and visualize the newly created model (where $S_1 = S$, $S_2 = I$, and $S_3 = R$).

sir_model = decoration(F_epi(sir_wiring_diagram));
Graph(sir_model)

Time domain

tmax = 40.0
tspan = (0.0,tmax);
(0.0, 40.0)

Initial conditions

u0 = [990,10,0]; # S,I,R
3-element Array{Int64,1}:
 990
  10
   0

Parameter values

p = [0.05*10.0/sum(u0),0.25]; # β*c/N,γ
2-element Array{Float64,1}:
 0.0005
 0.25

Random number seed

We set a random number seed for reproducibility.

Random.seed!(1234);
Random.MersenneTwister(UInt32[0x000004d2], Random.DSFMT.DSFMT_state(Int32[-
1393240018, 1073611148, 45497681, 1072875908, 436273599, 1073674613, -20437
16458, 1073445557, -254908435, 1072827086  …  -599655111, 1073144102, 36765
5457, 1072985259, -1278750689, 1018350124, -597141475, 249849711, 382, 0]),
 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 
0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 
0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0
0000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000
0000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000
0000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000
0000000000000000000000  …  0x00000000000000000000000000000000, 0x0000000000
0000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000
0000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000
0000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000
0000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000000
0000000000], 1002, 0)

Generating and running models

As ODEs

prob_ode = ODEProblem(sir_model,u0,tspan,p)
sol_ode = solve(prob_ode, Tsit5());
plot(sol_ode)

As SDEs

prob_sde = SDEProblem(sir_model,u0,tspan,p)
sol_sde = solve(prob_sde,LambaEM());
ERROR: MethodError: no method matching solve(::Tuple{DiffEqBase.SDEProblem{Array{Int64,1},Tuple{Float64,Float64},true,Array{Float64,1},Nothing,DiffEqBase.SDEFunction{true,Petri.var"#f#3"{Array{Int64,1},Array{Tuple{Dict{Int64,Int64},Dict{Int64,Int64}},1},Dict{Any,Any}},Petri.var"#noise#12"{Array{Tuple{Dict{Int64,Int64},Dict{Int64,Int64}},1},Dict{Any,Any},Dict{Int64,Int64},Dict{Int64,Int64}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Petri.var"#noise#12"{Array{Tuple{Dict{Int64,Int64},Dict{Int64,Int64}},1},Dict{Any,Any},Dict{Int64,Int64},Dict{Int64,Int64}},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SparseArrays.SparseMatrixCSC{Float64,Int64}},DiffEqBase.CallbackSet{Tuple{DiffEqBase.ContinuousCallback{Petri.var"#4#6"{Int64},Petri.var"#5#7"{Int64},Petri.var"#5#7"{Int64},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},DiffEqBase.ContinuousCallback{Petri.var"#4#6"{Int64},Petri.var"#5#7"{Int64},Petri.var"#5#7"{Int64},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},DiffEqBase.ContinuousCallback{Petri.var"#4#6"{Int64},Petri.var"#5#7"{Int64},Petri.var"#5#7"{Int64},typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}},Tuple{}}}, ::StochasticDiffEq.LambaEM{true})
Closest candidates are:
  solve(!Matched::DiffEqBase.EnsembleProblem, ::Any...; kwargs...) at /home/simon/.julia/packages/DiffEqBase/JNliX/src/solve.jl:125
  solve(!Matched::DiffEqBase.AbstractNoiseProblem, ::Any...; kwargs...) at /home/simon/.julia/packages/DiffEqBase/JNliX/src/solve.jl:133
  solve(!Matched::DiffEqBase.AbstractJumpProblem, ::Any...; kwargs...) at /home/simon/.julia/packages/DiffEqBase/JNliX/src/solve.jl:137
  ...
plot(sol_sde)
ERROR: UndefVarError: sol_sde not defined

As jump process

prob_jump = JumpProblem(sir_model, u0, tspan, p)
sol_jump = solve(prob_jump,SSAStepper());
plot(sol_jump)

Appendix

Computer Information

Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, icelake-client)
Environment:
  JULIA_NUM_THREADS = 4

Package Information

Status `~/.julia/environments/v1.4/Project.toml`
[80f14c24-f653-4e6a-9b94-39d6b0f70001] AbstractMCMC 1.0.1
[537997a7-5e4e-5d89-9595-2241ea00577e] AbstractPlotting 0.12.3
[46ada45e-f475-11e8-01d0-f70cc89e6671] Agents 3.2.1
[4f99eebe-17bf-4e98-b6a1-2c4f205a959b] AlgebraicPetri 0.3.1
[f5f396d3-230c-5e07-80e6-9fadf06146cc] ApproxBayes 0.3.2
[c52e3926-4ff0-5f6e-af25-54175e0327b1] Atom 0.12.16
[fbb218c0-5317-5bc6-957e-2ee96dd4b1f0] BSON 0.2.6
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.5.0
[a134a8b2-14d6-55f6-9291-3336d3ab0209] BlackBoxOptim 0.5.0
[2d3116d5-4b8f-5680-861c-71f149790274] Bridge 0.11.3
[1aa9af3a-2424-508f-bb7e-0626de155470] BridgeDiffEq 0.1.0
[46d747a0-b9e1-11e9-14b5-615c73e45078] BridgeSDEInference 0.3.2
[336ed68f-0bac-5ca0-87d4-7b16caf5d00b] CSV 0.7.3
[49dc2e85-a5d0-5ad3-a950-438e2897f1b9] Calculus 0.5.1
[134e5e36-593f-5add-ad60-77f754baafbe] Catlab 0.7.1
[aaaa29a8-35af-508c-8bc3-b662a17a0fe5] Clustering 0.14.1
[2445eb08-9709-466a-b3fc-47e12bd697a2] DataDrivenDiffEq 0.3.1
[a93c6f00-e57d-5684-b7b6-d8193f3e46c0] DataFrames 0.21.4
[7806a523-6efd-50cb-b5f6-3fa6f1930dbb] DecisionTree 0.10.6
[bcd4f6db-9728-5f36-b5f7-82caef46ccdb] DelayDiffEq 5.24.1
[2b5f629d-d688-5b77-993f-72d75c75574e] DiffEqBase 6.40.7
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 2.16.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 4.3.0
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.13.3
[aae7a2af-3d4f-5e19-a356-7da93b79d9d0] DiffEqFlux 1.17.0
[c894b116-72e5-5b58-be3c-e6d8d4ac2b12] DiffEqJump 6.9.3
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.16.0
[41bf760c-e81c-5289-8e54-58b1f1f8abe2] DiffEqSensitivity 6.23.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.15.0
[b4f34e82-e78d-54a5-968a-f98e89d6e8f7] Distances 0.9.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.23.4
[634d3b9d-ee7a-5ddf-bec9-22491ea816e1] DrWatson 1.14.4
[f6006082-12f8-11e9-0c9c-0d5d367ab1e5] EvoTrees 0.4.9
[587475ba-b771-5e3f-ad9e-33799f191a9c] Flux 0.10.4
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.12
[38e38edf-8417-5370-95a0-9cbb8c7f171a] GLM 1.3.9
[28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71] GR 0.50.1
[891a1506-143c-57d2-908e-e1f8e92e6de9] GaussianProcesses 0.12.1
[ea4f424c-a589-11e8-07c0-fd5c91b9da4a] Gen 0.3.5
[523d8e89-b243-5607-941c-87d699ea6713] Gillespie 0.1.0
[e850a1a4-d859-11e8-3d54-a195e6d045d3] GpABC 0.1.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.21.2
[a98d9a8b-a2ab-59e6-89dd-64a1c18fca59] Interpolations 0.12.10
[c8e1da08-722c-5040-9ed9-7db0dc04731e] IterTools 1.3.0
[4076af6c-e467-56ae-b986-b466b2749572] JuMP 0.21.3
[e5e0dc1b-0480-54bc-9374-aad01c23163d] Juno 0.8.2
[b1bec4e5-fd48-53fe-b0cb-9723c09d164b] LIBSVM 0.4.0
[b964fa9f-0449-5b57-a5c2-d3ea65f4040f] LaTeXStrings 1.1.0
[2ee39098-c373-598a-b85f-a56591580800] LabelledArrays 1.3.0
[23fbe1c1-3f47-55db-b15f-69d7ec21a316] Latexify 0.13.5
[7acf609c-83a4-11e9-1ffb-b912bcd3b04a] LightGBM 0.3.1
[093fc24a-ae57-5d10-9952-331d41423f4d] LightGraphs 1.3.3
[30fc2ffe-d236-52d8-8643-a9d8f7c094a7] LossFunctions 0.6.2
[c7f686f2-ff18-58e9-bc7b-31028e88f75d] MCMCChains 4.0.1
[add582a8-e3ab-11e8-2d5e-e98b27df1bc7] MLJ 0.12.0
[094fc8d1-fd35-5302-93ea-dabda2abf845] MLJFlux 0.1.2
[6ee0df7b-362f-4a72-a706-9e79364fb692] MLJLinearModels 0.5.0
[d491faf4-2d78-11e9-2867-c94bc002c0b7] MLJModels 0.11.0
[1914dd2f-81c6-5fcd-8719-6d5c9610ff09] MacroTools 0.5.5
[5424a776-8be3-5c5b-a13f-3551f69ba0e6] Mamba 0.12.4
[ff71e718-51f3-5ec2-a782-8ffcbfa3c316] MixedModels 3.0.0-DEV
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 3.13.0
[6f286f6a-111f-5878-ab1e-185364afe411] MultivariateStats 0.7.0
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.6.0
[9bbee03b-0db5-5f46-924f-b5c9c21b8c60] NaiveBayes 0.4.0
[b8a86587-4115-5ab1-83bc-aa920d37bbce] NearestNeighbors 0.4.6
[41ceaf6f-1696-4a54-9b49-2e7a9ec3782e] NestedSamplers 0.4.0
[47be7bcc-f1a6-5447-8b36-7eeeff7534fd] ORCA 0.4.0
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.21.0
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.41.0
[42b8e9d4-006b-409a-8472-7f34b3fb58af] ParallelKMeans 0.1.8
[4259d249-1051-49fa-8328-3f8ab9391c33] Petri 1.1.0
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 1.5.4
[c3e4b0f8-55cb-11ea-2926-15256bba5781] Pluto 0.10.6
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.9.0
[1a8c2f83-1ff3-5112-b086-8aa67b057ba1] Query 0.12.3-DEV
[6f49c342-dc21-5d91-9882-a32aef131414] RCall 0.13.7
[e6cf234a-135c-5ec9-84dd-332b85af5143] RandomNumbers 1.4.0
[c5292f4c-5179-55e1-98c5-05642aab7184] ResumableFunctions 0.5.1
[37e2e3b7-166d-5795-8a7a-e32c996b4267] ReverseDiff 1.2.0
[3646fa90-6ef7-5e7e-9f22-8aca16db6324] ScikitLearn 0.6.2
[f5ac2a72-33c7-5caf-b863-f02fefdcf428] SemanticModels 0.3.0
[428bdadb-6287-5aa5-874b-9969638295fd] SimJulia 0.8.0
[05bca326-078c-5bf0-a5bf-ce7c7982d7fd] SimpleDiffEq 1.1.0
[276daf66-3868-5448-9aa4-cd146d93841b] SpecialFunctions 0.10.3
[5a560754-308a-11ea-3701-ef72685e98d6] Splines2 0.1.0
[2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91] StatsBase 0.33.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.14.6
[789caeaf-c7a9-5a7d-9973-96adeb23e2a0] StochasticDiffEq 6.24.0
[92b13dbe-c966-51a2-8445-caca9f8a7d42] TaylorIntegration 0.8.3
[9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c] Tracker 0.2.8
[fce5fe82-541a-59a6-adf8-730c64b5f9a0] Turing 0.13.0
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 1.3.0
[276b4fcb-3e11-5398-bf8b-a0c2d153d008] WGLMakie 0.2.5
[29a6e085-ba6d-5f35-a997-948ac2efa89a] Wavelets 0.9.2
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.10.2
[009559a3-9522-5dbb-924b-0b6ed2b22bb9] XGBoost 1.1.1